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Abstract 

The changes of the solid Earth in south central Alaska in response to two major glacial 
fluctuations on different temporal and spatial scales have been estimated and we evaluated their 
influence on the stress state and ongoing tectonic deformation of the region. During the recent 
(1993-1995) Bering Glacier surge, a large transfer of ice from the Bagley Ice Field to the Bering 
Glacier terminus region occurred. We estimated the elastic displacement of the solid Earth due 
to ice mass redistribution from Global Positioning System (GPS) measurements at sites near the 
surging glacier. We can account for these displacements by transfer of an ice volume of about 14 
km from the surge reservoir area to the terminus region. We examined the background 
seismicity {Mi > 2.5) before, during, and after the surge. We found that the occurrence of small 
earthquakes (M L < 4.0) in the surge reservoir region increased during the surge time interval 
possibly in response to a decrease in ice mass. This suggests that a small decrease in the vertical 
stress, <j 3 , could be enough to modulate the occurrence of small, shallow earthquakes in this 
dominantly thrust fault setting. During this century the southern Alaska coastal glaciers have 
been undergoing an overall decrease in volume. Based on our compilation of changes in the 
extent and thickness of the coastal glaciers between the Malaspina and Bering, we calculated 
surface displacements due to the Earth’s viscoelastic response to annual thinning and to the 
cumulative retreat over the last 100 years. The uplift of the region due to an average annual 
thinning rate of 1-6 m/yr in the ablation region is 1-12 mm/yr. For our reference model with a 
viscosity of 5 x 10 19 Pa s for depths between « 40 and 200 km the total viscoelastic response due 
to the retreat over the last century may be as much as a couple of meters within the coastal 
ablation zone near Icy Bay. The maximum decrease in oy between 0 and 10 km was ssl.O MPa, 
which is significant in relation to the stress drops in recent earthquakes («2 to 10 MPa) but 
small in relation to the estimated tectonic stress magnitude. Therefore the occurrence of an 
earthquake such as the St. Elias (1979, Ms = 7.2) may have been advanced in time; however, 
most of the ongoing stress accumulation would be primarily due to tectonic forces. 
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1. Introduction 

The eastern Chugach Mountain range of southern 
Alaska is covered with a continuous series of con- 
nected glaciers (Figure 1) [Field, 1975]. Although 
individual glacier fluctuations axe variable and asyn- 
chronous, there has been a gross regional pattern of 
glacier retreat in southern Alaska this century [Afeier, 
1984; Porter , 1989; Molnia and Post , 1995]. By com- 
paring the predicted elastic response of the Earth to 
geodetic observations, the change in ice sheet mass 
can be estimated [Hager, 1991; Hager et al., 1991; 
Sauber et al , 1995; James and Ivins , 1995, 1998; 
Wahr et al., 1995]. Time-dependent deformation 
due to the viscoelastic response of the Earth to un- 
loading during the last 100 years is also likely to 
be significant. Here we report new constraints on 
the retreat of coastal glaciers this century between 
the Malaspina and Bering Glaciers. We estimated 
the magnitude of viscoelastic displacements associ- 
ated with this ice mass unloading and compared it 
with measured geodetic and longer-term deformation 
rates. We contrasted the importance of tectonic and 
glacial rebound in explaining the observed deforma- 
tion rate near Icy Bay. 

The recession of the Bering Glacier has been inter- 
rupted by at least six surges this century [e.g., Mol- 
nia and Post , 1995; Muller and Fleisher , 1995]. These 
surges involve periodic rapid movement of large quan- 
tities of ice within a glacier alternating with much 
longer periods of near stagnation or retreat [Afeier 
and Post , 1969; Molnia , 1993; Budd and Mclnnes , 
1974]. When a surge removes ice from the upper 
reaches of the glacier, its surface lowers by tens or 
hundreds of meters as ice is transported down glacier, 
where the ice thickens. Sometimes this is accompa- 
nied by an advance of the glacier terminus. The ice 
mass changes result in uplift of the solid Earth near 
the unloading (surge reservoir region) and subsidence 
beneath and near the receiving area [Cohen, 1993; 
Sauber et al., 1995]. In this paper we employ precise 
geodetic measurements made with the Global Posi- 
tioning System (GPS) at sites adjacent to the Bagley 
Ice Field and near the Bering Glacier (Figure 1) to 
supplement glaciological data to constrain ice mass 
redistribution, to estimate the total ice mass trans- 
fer and to explore the implications of our results for 
understanding the surge cycle of the Bering Glacier. 

The predicted stress changes associated with the 
Bering Glacier surge and glacier retreat in the last 100 
years are small in comparison with the tectonic stress 


levels estimated from borehole breakout data at a 
comparable depth. Based on earlier (water) reservoir- 
induced earthquake studies, however, we postulated 
that the small stress changes associated with glacial 
fluctuations this century could exert discernible con- 
trol on the occurrence of earthquakes. Seismicity 
(M l > 2.5) bracketing the time of the Bering Glacier 
surge and the location of moderate to large earth- 
quakes (M l > 4.0) between 1973 and 1997 were ex- 
amined to evaluate whether glacial fluctuations had 
discernible influence on earthquake occurrence. 

2. Representation of the Crust and 
Upper Mantle Rheological Properties 
for Estimating Surface Deformation 

The specific response of the Earth to a change in 
surface load across different spatial scales (10°-10 4 
km) and on a variable time scale (10°-10 4 years) 
depends on the rheological structure of the crust 
and mantle. The surface displacements associated 
with the recession of continental scale ice sheets [e.g., 
Peltier and Andrews, 1976], other Alaska glaciers, 
[e.g., Clark, 1977], and lake loads [e.g., Bills et al . , 
1994] have been used both to probe Earth rheology 
and to provide constraints on the unloading history. 

As was summarized by Kirby [1985], the thick and 
mechanically heterogeneous continental crust presum- 
ably plays an active role in determining the style of 
near-surface deformation. In addition, mobile aque- 
ous fluids are thought to play a major role, via pore 
pressure effects on brittle materials, in controlling 
rock strength of the shallow and midcontinental crust. 
The gradual transition from localized deformation 
along faults or fractures to distributed strain within 
the crust is thought to occur above a certain tem- 
perature, but the ratio of the least compressive stress 
(cr 3 ) to the differential stress (a = or i — cr 3 ) also is 
hypothesized to play an important role. 

In this study we calculated the response of the 
Earth to two glacial fluctuations. The Bering Glacier 
surge caused large localized surface changes over the 
time interval of 2 years. Since the spatial scale of the 
ice thickness changes is less than 30 km 2 , the crust 
will deform primarily elastically and little viscous de- 
formation is assumed to occur. The representation 
of Earth rheology that we used to calculate surface 
displacement was a layered elastic model, and the ge- 
ometrical complexity of the load changes were repre- 
sented by disks with a diameter of 5 km. 

For ice thickness changes over the wider region as- 
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sociated with glacial retreat in the last 100 years, the 
viscoelastic response of the upper mantle, and possi- 
bly the lower crust, needed to be calculated. Whether 
the upper mantle is best represented as a Newtonian 
or non-Newtonian fluid depends on the dominant de- 
formation mechanism [ Karato and Wu, 1993], Here 
the Earth’s response to a load-induced perturbation 
to ongoing tectonic processes is represented, and we 
assumed a linear viscous rheology. The lower crust 
and the upper mantle have effective viscosities that 
have been estimated to range from 10 18 to 10 22 Pa s 
with Maxwell times (ratio of viscosity to shear mod- 
ulus) of months to thousands of years. 

Lateral asthenospheric viscosity variations would 
presumably play an important role in calculations of 
deformation rates [Kaufman and Wu, 1998]. Thus we 
used a finite element model with a subducting slab 
[after Cohen, 1996] to represent the complex rheolog- 
ical structure at this plate boundary. Additionally, 
different glaciers between the Malaspina and Bering 
have variable retreat profiles. In this study, the vis- 
coelastic response of the Earth was calculated for a 
general two-dimensional retreat profile. The range of 
viscosity values tried was derived from work on post- 
seismic deformation in Alaska and lake and glacial 
unloading studies. 

3. Global Positioning System 

During June of 1993 and 1995, GPS measurements 
were made at the sites shown in Figure 1 for 1-12 
days. In 1997, only the coastal sites between Icy Bay 
and the Bering Glacier were observed (Table 1). Most 
daily observing periods were greater than 8 hours. 
GAMIT software [King and Bock, 1997] and GPS 
phase observations were used to estimate station co- 
ordinates, orbit, Earth orientation, and atmospheric 
parameters each day as described by Feigl et af.[1993]. 
We then used the GLOBK software [Herring, 1997] to 
estimate station coordinates and a velocity over some 
time period by combining these estimates and their 
covariance matrices with those from a similar analysis 
performed at the Scripps Orbital and Permanent Ar- 
ray Center (SOPAC) [Fang and Bock, 1995] using ob- 
servations from 30-60 stations of the global tracking 
network coordinated by the International GPS Ser- 
vice (IGS) for Geodynamics. The reference frame was 
defined by minimizing the adjustments in velocities of 
12 IGS stations, including Fairbanks, from their val- 
ues in the North American fixed International Ter- 
restrial Reference Frame (ITRF96) [Boucher et ai, 


1996]. 

The horizontal velocities of the stations from our 
study region are given in Table 1 in a North Amer- 
ican frame, obtained by rotating from the no-net- 
rotation frame of ITRF96 to North America using the 
NUVEL-1A global plate model [DeMets et al., 1994], 
This reference frame is most useful for comparing a 
tectonic model of ongoing deformation to our geodetic 
observations. The vertical velocities are given, how- 
ever, relative to Cape Yakataga to provide a regional 
reference frame. Of the sites given in Figure 1, Cape 
Yakataga is furthest from ice fluctuations, and the 
predicted tectonic uplift is small (i5 mm/yr or less). 
In a study of earlier GPS results which included data 
for 1993 and 1995 [Sauber et ai, 1997], we suggested 
that the daily scatter in horizontal position estimates 
from an individual survey fails by a factor of about 2 
to account for the errors with correlation times of sev- 
eral years. Based on an analysis of the daily position 
of globally distributed continuous GPS data, Mao et 
ai [1999] suggest that the formal error in the vertical 
component may be underestimated by a factor of 5 or 
greater. In this study the daily vertical repeatabilities 
show greater scatter than the horizontal components; 
so the formal errors have been scaled by a factor of 3. 

4. Tectonic Strain Accumulation 

In our study region, tectonic strain accumulation 
is due primarily to subduction of the Pacific plate and 
collision of the Yakutat terrane with interior Alaska 
[e.g., Plafker et ai, 1994]. In this geologically complex 
region between the transcurrent Fairweather fault and 
the Alaska- Aleutian subduction zone, recent crusted 
shortening and strike-slip faulting occurs offshore in 
the Gulf of Alaska (1987-1988, M s = 6.9, 7.6, 7.6) and 
onshore in the Chugach-St. Elias Mountains (1979, 
Ms — 7.2). Prior great earthquakes in the region 
occurred in 1899 ( Mw = 8.1, Yakataga; Mw = 8.1, 
Yakutat Bay) [Thatcher and Plafker, 1977, unpub- 
lished manuscript, 1977]. 

The tectonic process assumed to exert the greatest 
influence on the geodetic observations reported in this 
study is deformation associated with a locked plate 
interface at shallow depths (<40 km) [Savage and 
Lisowksi, 1988]. The horizontal rate of deformation at 
stations located more than 20 km from major glacial 
fluctuations are consistent with the deformation rate 
predicted from elastic dislocation models of a locked 
main thrust zone [Sauber et ai, 1997], Of the stations 
given in Figure 1 and Table 1, only the Yakataga sta- 
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tions (VYAK, YAKU, FURR) and, to a lesser extent, 
TIME were used to constrain interseismic strain ac- 
cumulation models. Relative to stable North Amer- 
ica the horizontal rates of tectonic deformation across 
the region given in Figure 1 axe predicted to be 20-40 
mm/yr, and the estimated tectonic uplift rates range 
from approximately zero near the Gulf of Alaska coast 
up to 12 mm/yr further inland. 

5. Crustal Deformation in the Region 
Near the Surge of the Bering Glacier 

The GPS-derived displacement between 1993 and 
1995 reflects crustal deformation due to tectonic and 
nontectonic forces. In the previous section we dis- 
cussed the predicted short-term tectonic strain, and in 
this section we discuss the crustal deformation due to 
large changes in ice thickness as a result of the 1993- 
1995 glacier surge. The GPS measurements cannot be 
used to uniquely constrain both the location and mag- 
nitude of ice thickness changes. Therefore glaciologi- 
cal data were used to estimate the general region that 
underwent ice thickness changes, as well as to provide 
some constraints on the relative magnitude of these 
changes (Plate 1). The GPS results were then used 
to test alternate ice transfer models suggested from 
the glaciological data. Our initial ice change model 
prompted us to make an aircraft flight over the surge 
reservoir region (B. Molnia, August 1999). We were 
able to identify trimlines, especially on south and east 
facing slopes (for example, near station Isle, Plate 2), 
associated with ice thinning attributed primarily to 
the surge. 

5.1. Glaciological Constraints on Transfer of 
Ice Mass During a Surge 

The results of an extensive effort to study the 1993- 
1995 Bering Glacier surge have provided some con- 
straints on its timing, spatial extent, and ice thick- 
ness changes [Lingle et al., 1993; Molnia, 1993; Mol- 
nia et al., 1994; Molnia and Post, 1995; Roush, 1996; 
Herzfeld and Mayer, 1997; Fatland, 1998]. The surge 
seems to have originated south of the equilibrium line 
in the spring of 1993. Rapid ice movement down- 
glacier into the piedmont lobe and up-glacier into the 
Bagley Ice Field followed [Lingle et al., 1993; Fatland, 
1998], By late in the summer of 1993 the terminus 
began to advance. Ice transfer to the receiving area 
resulted in terminus advance of about 5 km along its 
30-km-wide front [Krimmel, 1994], and parts of the 
terminus advanced approximately 9 km [Molnia et al., 


1994]. 

Fatland and Lingle [1998] and Fatland [1998] used 
Cband synthetic aperture radar (SAR) interferome- 
try to estimate surface ice velocities on the Bagley 
Ice Field prior to and during the surge. Their studies 
documented regions of fast moving ice in the eastern 
and western Bagley Ice Field that extended up to ele- 
vations of about 1500 m. Additionally, aircraft flights 
over the region provided some constraints. Stranded 
snow on the valley wall 25—100 m above the drawdown 
of the lower Bagley Ice Field suggested extensive low- 
ering (B. Molnia, field observations, 1993, 1994). In 
August 1999 we took photographs and videotaped the 
ice margin during a fixed wing aircraft flight over the 
portion of the glacier involved in the 1993-1995 surge 
(B. Molnia, 1999). For example, the changes shown 
by Plate 2 reflect the cumulative thinning due pri- 
marily to the surge as well as annual thinning in the 
ablation zone. 

5.2. Model of Ice Thickness Changes 

On the basis of glaciological field observations, we 
identified the general region that underwent thinning. 
We specified vertical ice lowering in the surge reser- 
voir over a broad region between elevations of «900 
m on the upper reaches of the Bering Glacier and 
£31500 m on the Bagley Ice Field. On the basis of ice 
velocities during the surge we created a general rel- 
ative unloading model. Specifically, ERSl synthetic 
aperture radar data from the winter of 1994 were used 
by Fatland [1998] to estimate surface, horizontal ice 
velocities of 0.3 m/d in the eastern Bagley Ice Field 
at elevations of £31500 m and up to 4.5 m/d in the 
Bagley Ice Field near the top of the Bering Glacier 
at £31220 m. The region with fast moving ice during 
the surge was at elevations below the equilibrium line 
altitude (ELA). We assumed the region of greatest 
extension and thinning is associated with the highest 
surge ice velocities in the reservoir region, and we ta- 
pered the thickness change to zero at the limit of fast 
moving ice (Figure 19 of Fatland [1998]). 

This unloading model included both surge-related 
ice thinning and annual thinning associated with re- 
treat. Near the Tana Glacier, retreat of the glacier is 
the primary source of ice thinning. For other parts 
of the ice thinning region given in Plate 1 the annual 
thinning is as much as an order of magnitude smaller 
than the surge change. The value assumed for the ini- 
tial unloading model over the 1993-1995 time interval 
shown in Plate 1 is given from top left to bottom right 
in meters: -15, -15, -10, -30, -45, -45, -30, -30, -15, -5, 
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-5, -5, -35, -35, -10, -20, -15, -5, -15, -15, -10, -10, -5. 

The magnitude of ice thickening during the surge 
was greatest in the region in which the glacier ad- 
vanced (Plate 1). Additionally, thickening of the 
Bering Glacier piedmont lobe was estimated to be 
40-150 m with mobile bulges exceeding 200 m in 
thickness [tfous/i, 1996; B. Molnia, field observations, 
1993-1996]. We used ERS1 data from the Bering 
Glacier piedmont lobe to identify regions of ice ad- 
vance and thickening (presurge position definition is 
from images on June 16, 1992, and April 20, 1993; 
postsurge, from a September 22, 1995, image). 

Since the water load associated with Vitus Lake 
(northwest of the station DON) was replaced by a 
thicker ice load, we did place some disks in this region, 
but we made them thinner than the disks due north 
and northeast of the station DON. Unfortunately, we 
had just one geodetic station in this region, and we 
are unsure of the reliability of the displacement for 
estimating ice load changes; other processes such as 
sediment loading offshore due to the high flux of sed- 
iments associated with the surge could have been im- 
portant as well. 

An initial relative model of ice loading is given from 
top left to bottom right in meters: 10, 10, 10, 10, 10, 
10, 10, 10, 10, 10, 30, 30, 45, 40, 50, 50, 50. 

We made the simplifying assumption that the ice 
volume removed (L“) during 1993-1995 surge is ap- 
proximately equal to the ice volume added (L + ) to the 
Bering Glacier piedmont lobe. That is, J2iL\ L = 
Lj+i where p ice xH, pi ce equals the den- 
sity of ice, H corresponds to the ice thickness due to 
the surge (red only; black in the plate is attributed 
primarily to annual retreat of the Tana) change over 
a given time period, i equals the number of disk loads 
to represent the change in ice thickness in the Bagley 
Ice Field and upper reaches of the Bering Glacier, and 
j equals the number in the Bering Glacier piedmont 
lobe (Plate 1, blue values). This assumption is sup- 
ported by observations of ice transfer during other 
surges in Alaska and in the Pamirs of Asia [Dolgushin 
and Osipova , 1975; Kamb et al , 1985]. 

5.3. Elastic Displacements Caused by 
Redistribution of Ice Mass 

For an individual disk load (Li) of density p tce and 
radius a, the vertical (14) and horizontal (t/*) displace- 
ment as a function of distance from the center of a 
disk are given in terms of hypergeometric functions 
by equations 12-16 of Farrell [1972]. We assumed 


that the instantaneous response of the solid Earth to 
ice thickness changes over the small spatial aperture 
of the surge region was primarily elastic [Sau&er et 
al, 1995]. To represent some of the spatial complex- 
ity of a variable surface load, we chose to represent 
the change in glacial load by multiple disk loads each 
having a 5-km diameter. To approximate an equiva- 
lent rectangular load, the disks were multiplied by a 
geometrical scaling factor. The total displacement, u, 
v, at an individual geodetic station is the sum of the 
contributions from n disk loads Li , u «> Z)”=i v i- 

The vertical displacement of the solid Earth, u, is 
useful for estimating the magnitude (and sign) of ice 
thickness changes. The horizontal component is par- 
ticularly sensitive to the location (direction) of load 
changes. 

5.4. Comparison of Observed and Predicted 
Displacements 

Alternative models were tested against the 1993- 
1995 GPS displacement values by scaling the initial 
distribution of ice thickness change values above. We 
most closely matched the GPS results and remained 
consistent with the observed trimlines, with a scaling 
factor of 1.4. A set of unloading/loading disks that 
can account for the estimated station displacements is 
given in Plate 1 and the predicted horizontal and ver- 
tical displacements are shown in Plate 1 and Figure 
2. The maximum thinning is estimated to be about 
63 m, and the maximum ice thickening is 70 m. The 
predicted vertical depression ranges up to 102 mm in 
the region of increased ice thickness in the terminus 
region and the predicted uplift ranges up to 94 mm in 
the unloading region where the Bagley Ice Field flows 
into the Bering Glacier. 

In Table 2 we compared the displacements pre- 
dicted by our best fitting forward model to the geode- 
tic observations. The observed displacement has 
had the estimated tectonic displacement removed [see 
Sauber et al., 1997]. Of the sites used in this surge 
study, just DON had a third set of observations which 
enabled us to independently estimate the tectonic 
component of the displacement based on the 1995- 
1997 results. For ANCX, ISLE, and TIME the tec- 
tonic deformation rate had to be estimated on the 
basis of tectonic modeling of geodetic results from sta- 
tions farther (>30 km) from the surge, and some addi- 
tional uncertainty (±5 mm/yr) should be attributed 
to this approximation. 

At the three sites nearest the surge (DON, ANCX, 
and ISLE) our model predictions minus the observed 
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(corrected) vertical displacement over the 1993-1995 
time interval are less than the formal GPS uncertainty 
in the measured values (wlO mm). The horizontal 
component is very sensitive to the direction of the 
ice mass change (i.e., changes in azimuth to unload- 
ing/loading will result in significant changes in magni- 
tude and sign of the displacement). Because of limited 
glaciological and geodetic observations, we could use 
only a moderate resolution grid; this made it difficult 
to model the horizontal geodetic data in particular. 
We do not understand the large discrepancy between 
the observed and predicted north value for the station 
DON; the geodetic results seem to imply an ice mass 
change south of the stations. 

The glaciological data (aerial photographs, field 
studies, and interferometric SAR study of ice veloci- 
ties) were the primary constraints on the redistribu- 
tion of ice mass, and we axe confident of the general 
region that underwent thinning and thickening. The 
GPS data provided the most useful information near 
the stations ISLE and ANCX; here the ice thickness 
change uncertainty may be ±10-20 m (see the sensi- 
tivity analysis summarized by Sauber et al. [1997]). 
In other areas the uncertainties are larger and the 
magnitude and distribution of change were estimated 
primarily by glaciological studies. 

5.5. Postsurge Changes to Bering Glacier 
and Bagley Ice Field 

Between the cessation of the surge in 1995 and Au- 
gust 1999, the region that thickened and advanced 
during the surge has undergone thinning, and the 
terminus has retreated 0.1 to 1 km per year. The 
Bering Glacier piedmont lobe has a smooth surface 
once again. During the postsurge time period, most 
of the area below the ELA has experienced signifi- 
cant thinning. Although the surge initiation region 
may begin to build up ice thickness in the future ow- 
ing to transfer of material from upglacier, it is not 
evident at this point. 

6. Crustal Deformation Due to Glacial 
Recession This Century 

Nearly all of the glaciers shown in Figure 1 have 
receded during this century [e.g., Meier, 1984], and 
many of these glaciers are still undergoing retreat. To 
calculate the predicted crustal deformation over the 
time span of our geodetic observations, we needed to 
estimate the Earth response to recent changes as well 
as the cumulative retreat over the last 100 years. 


Here we report the available constraints on glacial 
retreat in the last 100 years. We calculated the sur- 
face deformation rate across the region due to an av- 
erage annual ice thinning rate assuming a simple lay- 
ered elastic Earth model. A two-dimensional finite 
element model was used to calculate the predicted 
displacement rate due to retreat. Because of our un- 
certain knowledge of the Earth’s response on the time 
scale of tens of years, we explored simple variations 
in Earth rheology (i.e., especially the asthenospheric 
viscosity) assuming a theoretical retreat profile. 

6.1. Glaciological Constraints 

Some of the data we used for a rough characteri- 
zation of ice thickness changes axe given in Table 3; 
they include the elevation of moraine crests and trim- 
lines in comparison to recent glacier surfaces. The Icy 
Bay glaciers have retreated ss35 km from their turn 
of the century terminal position at the mouth of Icy 
Bay [Plafker and Miller, 1958; Molnia, 1977; Porter, 
1989], and radiocarbon dating of lateral moraines sug- 
gest thinning of «300 m over this same time period 
(G. Plafker, field observations, 1963, 1969, 1982). The 
position of the Malaspina Glacier terminus is essen- 
tially stagnant. On the basis of trimline heights in the 
Samovar Hills and Chaix Hills, however, ss 140-180 m 
of thinning of the inner margin of the large piedmont 
lobe has occurred, probably during this century (G. 
Plafker, field observations, 1963, 1969, 1982; B. Mol- 
nia, field observations, 1974, 1989-1991, 1998). Thin- 
ning of the Bagley Ice Field and Tana Glacier during 
this century, but prior to the 1993—1995 surge, has 
been estimated to be 27 to 90 m at an average ele- 
vation of «1500 m [Miller, 1957; B. Molnia, field ob- 
servations, 1974-1993; J. Sauber, field observations, 
1993]. Between the 1967 surge and the onset of the 
most recent surge in 1993, the Bering Glacier termi- 
nus receded as much as 10 km and thinned as much 
as ss 180 m [Molnia, 1993; Molnia and Post, 1995]. 
The elevation changes along the centerline of Bering 
Glacier have been determined directly from a com- 
parison of the 1972 and 1990-1992 elevations (Figure 
lb of Molnia and Post [1995]). For the other glaciers 
in this region this level of detail is not available. 

Where we do not have a direct estimate of ice thin- 
ning, we used the change in the position of the glacier 
terminus to determine a rough value. On the basis of 
a study of 15 mountain valley glaciers, the ratio of the 
thickness change averaged over the full length to the 
change in terminus position has been characterized 
by a profile shape factor, /, estimated at 0. 1-0.4, with 
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an average value of ss0.3 [Schwitter and Raymond, 
1993]. It should be noted, however, that Sapiano et 
al. [1998] examined the elevation, volume, and termi- 
nus changes of nine glaciers in Alaska and Washing- 
ton and found no simple relationship between volume 
change and terminus retreat. However, they found 
the elevation changes in most cases were largest near 
the terminus and decreased upglacier more rapidly 
than a linear variation with distance. 

Field [1975] estimated the transition between the 
ablation and accumulation zone from the annual firn 
limit, a proxy for the EL A. He found that the aver- 
age was generally around 1000 m (3280 feet) near the 
coast and the limit rises to about 1500 m (4920 feet) 
inland. We did not place disks representing unloading 
above these elevations. 

6.2. Elastic Response to Average Yearly 
Thinning Rates 

As was discussed in the preceding section, to repre- 
sent some of the spatial complexity of a variable sur- 
face load, we chose to represent the change in glacial 
load by multiple disk loads. Since detailed informa- 
tion on retreat was not available for many parts of 
this region, we used a disk diameter of 10 km to rep- 
resent average changes at this spatial scale (Figure 3), 
and we used Farrell [1972] to calculate the elastic dis- 
placement assuming a spherical Earth model. Also, 
we had to assume a temporal average over the time 
period of the constraints given in Table 3. 

The calculated elastic uplift rate caused by the av- 
erage yearly load reductions are given in Figure 3 for 
the horizontal component and in Figure 4 for the ver- 
tical component. The distributed elastic response to 
changing ice loads is up to 12 mm/yr of uplift and 2 
mm/yr for the horizontal components. 

We note that the differences between the horizon- 
tal (north and east) and vertical velocities of the five 
coastal stations (VYAK, FURR, YAKU, AMBR, and 
DON (1995-1997)) and the weighted mean of the indi- 
vidual components for all sites are less than 2cr except 
for the north and vertical component at the station 
AMBR (near Icy Bay). The uplift of AMBR (Icy Bay 
region) relative to VYAK (Yakataga region) is about 
12 mm/yr (see Table 1), and the horizontal displace- 
ment rate relative to Yakataga is 6 mm/yr (mostly 
south). As can be seen in Figures 3 and 4, the pre- 
dicted annual displacement due to retreat near the 
station AMBR is vertical uplift (about 3 mm/yr) and 
southward displacement (about 1 mm/yr). 


As is evident in Figure 3 the orientation and mag- 
nitude of surface displacement rate are complex on 
the local and regional spatial scale. This illustrates 
the importance of good regional glaciological data to 
model displacement rates obtained from geodetic and 
tide gauge measurements. These calculations do not 
account for the viscoelastic response of the Earth to 
glacial unloading since early this century. 

6.3. Viscoelastic Displacement Associated 
With Retreat 

We used a two-dimensional plane strain finite ele- 
ment method (TECTON [Melosh and Raefsky, 1981]) 
to calculate the viscoelastic response of the Earth to 
glacier retreat during the past 100 years. The finite 
element grid across this subduction zone plate bound- 
ary is modified from Cohen [1996] and includes a 
shallow dipping subducting slab and both an oceanic 
crust-mantle and a continental crust-mantle (Figure 
5). Since the details of crustal and upper mantle 
lithology and temperature gradient are not available 
across the study region, we relied on an estimate of 
the transitional depth from the epicenters of large 
earthquakes and the depth of background seismicity 
[e.g., Scholz, 1990]. In south central Alaska, upper 
plate background seismicity extends to less than 40 
km [ Page et al., 1991], and the downgoing Pacific 
plate begins to bend at «22 km [Plafker et al., 1994], 
For the oceanic plate we assumed a transitional depth 
of around 30 km [Sauber et al, 1993]. 

We used viscosities of 10 25 Pa s for the upper crust 
(0-38 km) and 10 21 Pa s for the upper mantle between 
210 and 500 km. Owing to the uncertainty in the vis- 
cosity of the lower crust and upper mantle (>38 km 
to 210 km in our nonunique model for the continental 
plate and between 30 and 210 km below the oceanic 
plate), we calculated the predicted displacements by 
assuming a range of viscosities (5 x 10 18 Pa s to 5 x 
10 21 Pas). This range of numerical values is represen- 
tative of those found from studies of postseismic and 
glacial rebound at the subduction boundary in Alaska 
[Savage and Plafker, 1991; Cohen, 1996; Zheng et al, 
1996] and Cascadia [e.g., Wang et al, 1994; James et 
al, 2000] and from lake unloading studies in a tecton- 
ically active region [e.g., Bills et al, 1994]. This does 
not include a short-term, downdip low- viscosity zone 
or a creep zone used to account for rapid postseismic 
slip [e.g., Wahr and Wtyss, 1980; Cohen, 1996]. 

The uplift rate due to glacial retreat during this 
century will be largest in the ablation zones of glaciers 
near the Gulf of Alaska coast. As is seen in Figure 1, 
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there is variation in the location of the glaciers relative 
to the coast. The Bering Glacier piedmont lobe has 
undergone advance during periodic surges followed by 
slow retreat of the ice front between surges, and the 
changes in the Bering Glacier lobe may have a small 
viscoelastic response. On the other hand, between 
the Malaspina Glacier and the glaciers to the west 
of Icy Bay, significant retreat and ice thinning (tens 
to hundreds of meters) have occurred, and we would 
expect a higher rebound in this area then, for exam- 
ple, near Cape Yakataga (Figure 1). Thus there is 
a suggestion that near Icy Bay the crustal deforma- 
tion in response to the retreating glaciers may have 
a discernible influence on incremental tectonic strain 
data. 

A general idealized longitudinal unloading profile 
for a retreating valley glacier is given in Figure 6 (af- 
ter Figure 1 of Schwitter and Raymond [1993]). This 
roughly corresponds to a north-south profile that in- 
cludes the site AMBR (Figure 3) and is near a site 
with carbonl4 dating of an overridden forest. The 
estimated ice thickness change near this site is about 
300 m (Table 3, site 1). For the region including the 
peak in Figure 6, spanned by a 10-km element in the 
finite element grid, a maximum of 300 m of ice thin- 
ning was assumed; the other elements were scaled on 
basis of the general profile from Schwitter and Ray- 
mond [1993]. For simplicity, we assumed that most 
glacier retreat and thinning occurred 100 years ago. 
We then calculated the predicted horizontal and verti- 
cal displacement predicted for shortly after the retreat 
until 100 years later. 

For illustrative purposes, in Figure 7 we show the 
range of predicted deformation rates for 1, 50, and 
100 years for the asthenospheric values given in Fig- 
ure 5. The maximum uplift of about 3 m is centered 
near the Gulf of Alaska coast and drops to zero at 
«200 km inland from the coast; there is some migra- 
tion of the maximum uplift away from the coast with 
time. Note that the station AMBR is located near 
the region of maximum rebound. In addition to the 
uncertainly associated with the GPS measurement 
results, the position of this site with respect to the 
actual ice thinning profile is approximate. 

The range of deformation rates due to the vis- 
coelastic response to retreat are compared in Table 4 
with the uplift and horizontal rate of AMBR relative 
to VYAK. For asthenospheric viscosities >5 x 10 2 ° 
Pas, very little time-dependent deformation over the 
10 years that span our geodetic observations at the 
station AMBR is predicted. In contrast, with an as- 


thenospheric viscosity of 5 x 10 18 Pa s the vertical 
uplift rate due to the viscoelastic response was esti- 
mated to be about 31 mm/yr. A value of 5 x 10 19 
Pa s is most consistent with the observed rate, but 
the results from AMBR are not well determined, and 
our model contains simplifying assumptions. An up- 
lift rate of >1 cm/yr probably cannot be accounted 
for by the annual expected thinning alone. An as- 
thenospheric viscosity of less than 10 20 Pa s may be 
necessary to account for the preliminary GPS results. 

6.4. Crustal Deformation From Late 
Pleistocene Deglaciation 

The last major Pleistocene deglaciation episode 
in Alaska is correlated with late Wisconsin fluctua- 
tions of the Laurentide Ice Sheet. Radiocarbon ages 
the show onset of glaciation at ^24,000 B.P. and 
deglaciation beginning at about 13,500 B.P. [ Molnia , 
1989; Hamilton , 1994]. The late Pleistocene ice un- 
loading model ICE-4G includes a simple representa- 
tion of deglaciation in southern Alaska and was used 
to obtain a rough estimate of the solid Earth dis- 
placement caused by viscous relaxation [ Peltier , 1993, 
1994; T. James, personal communication, 1998]. An 
estimate of the present-day uplift attributable to this 
deglaciation was obtained by examining vertical dis- 
placements since 1000 B.P. obtained from files of to- 
pography change computed by W. R. Peltier. These 
files are available through the National Geophysical 
Data Center, Boulder, Colorado (see also the discus- 
sion by James and Ivins [1998]). Along a profile per- 
pendicular to the coast near Cape Yakataga (Figure 
1), the uplift rate caused by late Pleistocene deglacia- 
tion was estimated to be 2 mm/yr. 

The asthenospheric viscosity of less than 10 20 Pa s 
suggested for our study region is lower than that de- 
rived for global postglacial rebound associated with 
continental scale glacial retreat in the late Pleistocene 
[e.g., Peltier , 1994]. It should be noted that the up- 
lift data used in the later study were primarily from 
nontectonic regions, and much of the late Pleistocene 
deglaciation occurred at interior plate regions. The 
ICE-4G results thus give an upper bound on the up- 
lift rate in southern Alaska due to late Pleistocene 
deglaciation. With a lower asthenospheric viscosity 
the predicted uplift due to late Pleistocene deglacia- 
tion would be even smaller, since most of the rebound 
would have occurred much earlier. 
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6.5. Coastal Uplift Rates Between Icy Bay 
and Bering Glacier 

Table 1 presents estimates of the horizontal and 
vertical rates of deformation from GPS observations 
made between 1993 and 1997, and Table 5 gives the 
uplift rates estimated from carbon 14 dating of ter- 
races uplifted over the last 3000-6000 years. Since 
this is a region of large earthquakes, significant elas- 
tic strain is accumulating which will eventually be 
released in earthquakes, and the short-term strain is 
thought to be primarily elastic. These earthquakes 
could cause slip on the plate interface and/or on faults 
within the overriding plate, causing crustal shortening 
and permanent uplift. Coseismic slip on faults within 
the overriding plate is suggested by sal m of probable 
coseismic uplift between Icy Bay and Yakataga dur- 
ing the great 1899 Yakutat Bay earthquake sequence 
[Parr and Martin, 1912; Thatcher and Plafker, 1977, 
unpublished manuscript, 1977]. Estimates of the up- 
lift rate over the last «6000 years suggest 7-15 mm/yr 
of uplift near Icy Bay but only 2-3 mm/yr near the 
Bering Glacier (Table 3, Plafker et al. [1981] and 
Molnia and Post [1995]). Thus more onshore crustal 
shortening may occur in the Icy Bay region. Also, as 
was discussed above, it may be that in the recently 
deglaciated region of Icy Bay there may be some ad- 
ditional short-term uplift due to rebound. 

7. Glacial Fluctuations and 
Earthquakes 

The ice mass changes due to the glacial fluctua- 
tions discussed above perturb the local and regional 
stress field. The spatial (Plate 1) and temporal (1- 
2 years) scales, as well as the magnitude, of load- 
ing and unloading associated with the Bering Glacier 
surge is similar to water reservoir impoundment or 
the removal of rock in a quarry. For this case, the fi- 
nite strength of the lithosphere will support the load, 
and the elastic, deviatoric stresses will decrease as a 
function of depth [e.g., Scholz, 1990; Cohen, 1993], 
The retreat of the coastal Alaska glaciers during this 
century has occurred over a broader spatial scale, and 
some flow in the asthenosphere is assumed to have oc- 
curred. Eventually, this flow will bring the Earth back 
into isostatic equilibrium. In the glacial retreat case, 
there will be time-dependent changes in surface and 
subsurface stress distribution even if the load does not 
continue to change. 

We evaluated the possible influence of glacial fluc- 
tuations on earthquake occurrence and surface fault- 


ing by presenting the predicted stress changes in the 
context of a simple Navier-Coulomb failure criterion. 
For the surge case we tested if there was a concur- 
rent seismicity change in the area of the surge reser- 
voir region where ice thinning occurred and/or in 
the surge receiving area where ice thickening and ad- 
vance occurred. Since retreat of the coastal Alaska 
glacier started approximately 100 years ago, we can 
not easily test for temporal variations in earthquake 
occurrence. We compared the seismicity in recently 
deglaciated regions to regions at greater distance from 
retreating glaciers. 

7.1. Magnitude and Orientation of Principal 
Stresses 

Earthquake focal mechanisms, offshore in situ bore- 
hole measurements, the geodetic estimate of incre- 
mental strain, and regional geological evidence [ Es- 
tabrook and Jacob, 1991; Plafker et al., 1994; Doser 
et al., 1997; Sauber et al., 1997] have been used to es- 
timate principal shortening and stress directions. For 
the coastal region between Icy Bay and Kayak Island, 
a horizontal north-south to northwest-southeast ori- 
entation is suggested for the maximum effective stress 
(cti), and a minimum effective stress (<r 3 ) that is verti- 
cal has been assumed for the dominantly thrust fault- 
ing environment (Figure 1). Strike-slip faulting as- 
sociated with vertical intermediate stress (o 2 ) condi- 
tions has been suggested near the Contact fault just 
north of the Bagley Ice Field [ Savage and Lisowski, 
1988]. 

Borehole failure observed in offshore wells between 
Icy Bay and Kayak Island have been reported by 
Hottman et al. [1979]. In general, breakouts have 
their long axis parallel to the minimum horizontal 
stress so that they can be used to map horizontal 
principal stress trajectories and identify the relative 
horizontal stress magnitudes [ Adams and Bell, 1991]. 
On the basis of the prevalence of thrust faulting in the 
area and independent information on the orientation 
of the maximum horizontal stress direction, Hottman 
et al. [1979] estimated the orientation and the approx- 
imate magnitudes of principal stresses. In their calcu- 
lations it was assumed that the three principal Earth 
stresses were oriented almost vertically and horizon- 
tally. The overburden stress, taken to be cr„, was 
estimated from well interval density logs. They cal- 
culated Earth stress gradients of <ti = 32 kPa/m, a 2 
- 27 kPa/m, and a v = <73 = 23 kPa/m, and pore 
pressure p a = 19.2 kPa/m from cores cut between 
2700 and 4000 m. In Figure 8, the estimated effective 
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stress for o\ and <t 3 were extrapolated to a depth of 
5 km. 

7.2. Shear Failure 

The Mohr-Coulomb criterion for brittle shear fail- 
ure in rock is described by 

r = T 0 + (1) 

where r is the shear stress necessary to induce failure 
on a fault plane, r 0 is the inherent shear strength of 
the fault, fx is the coefficient of friction of the fault sur- 
face, and <r n is the normal stress on the fault [ Jaeger , 
1969]. The locus of shear (r) and normal (cr n ) stress 
components on a suite of faults of various orienta- 
tions is a Mohr circle whose center is the average be- 
tween the maximum (<7i) and minimum (<73) principal 
stresses and whose radius is (<7i - cr 3 )/2. 

Laboratory studies indicate for intact rock samples 
that r 0 « 50 MPa and p. = 0.6 to 0.85 (e.g., see the 
summary by Johnston [1987]). For a region such as 
southern Alaska that is heavily faulted, we assume 
that faults with fault gouge exist at a variety of ori- 
entations. Additionally, when preexisting faults with 
fault gouge are present, p, may be as low as 0.2-0. 4, 
r Q approaches zero, and equation (1) becomes 

r = (0.4 - 0.85)er n . (2) 

For faults with low friction, in fact, slip on faults over 
a range of orientations would occur. 

7.3. Stress Drop in Recent Earthquakes 

A number of large earthquakes have occurred in the 
last 30 years within or near the study region given in 
Figure 1 (1987-1988, Gulf of Alaska, M s = 6.9, 7.6, 
7.6; 1979, St. Elias, Ms = 7.2; 1970, Pamplona zone 
Mw — 6.7). The static stress drop in the St. Elias 
and Gulf of Alaska earthquakes ranges from moderate 
to high («2 to 10 MPa) [Hwang and Kanamori , 1992; 
Estabrook et al., 1992; Sauber et ai , 1993; Doser et 
ai, 1997]. Although stress levels on individual faults 
are highly variable, these earthquakes suggest that 
much of the region given in Figure 1 is close to fail- 
ure. The last Ms >8.0 earthquakes occurred in 1899. 
In Figure 8 we present the stress drop in these earth- 
quakes relative to other stress changes. 

7.4. Surge Case 

Water loading, quarry unloading, and changes in 
glacier mass cause a static change in load and pore 


pressure equivalent to the water or rock thickness 
(see the summary by Scholz [1990]). A comparison 
of well-documented case histories of seismicity near 
water reservoirs suggests that primarily two types of 
induced seismicity are observed (see the summary by 
Simpson et ai [1988]). Changes in seismicity which 
follow rapidly after the filling of a reservoir are re- 
lated to changes in elastic stress or changes in pore 
pressure coupled to the elastic stress. Since the stress 
increase from the elastic load drops off rapidly with 
distance, seismicity in these cases is concentrated in 
the immediate vicinity of the reservoir, and earth- 
quake sizes tend to be small, since the stress gradients 
are high. Simpson et ai [1988] further observed that 
where there is a long delay between the filling of the 
reservoir and the start of increased seismicity, diffu- 
sion of pore pressure from the reservoir to hypocentral 
depths may play a dominant role. 

In the dominantly thrust earthquake environment 
of the study region, the direct effect of ice loading in 
the surge receiving region will be to increase cr 3 , which 
reduces the likelihood of earthquake faulting, whereas 
for the unloading region it will decrease a 3 and make 
failure more likely. At shallow crustal depths the im- 
mediate poroelastic effect could also be important. 
However, after the cessation of the surge, the Bering 
Glacier started to retreat once again; so we did not 
evaluate the possibility of later triggering due to pore 
pressure diffusion. 

We examined the background seismicity (Ml > 
2.5) before (1991.0-1993.4), during (1993.4-1995.8), 
and after (1995.8-1998.3) the Bering Glacier surge to 
see if a short-term change in either the rate or loca- 
tion of small to moderate earthquakes occurred. Al- 
though the general regional level of seismicity varies 
only slightly over the three time periods, there is an 
increased level of seismicity under the region of great- 
est thinning during the surge (Figure 9). An ice thick- 
ness change of 50 m corresponds to a change in the cr 3 
(vertical) of approximately 0.5 MPa. This is consis- 
tent with the idea that 03 (vertical) is decreased (Fig- 
ure 8) and failure in thrust type earthquakes would be 
more likely; the focal mechanism solutions, however, 
for these small earthquakes are unknown. 

7.5. Retreat Case 

We used the plane strain viscoelastic calculations 
described in an earlier section to examine the esti- 
mated stress changes as a function of time and depth. 
The calculated values represent average values over 
the elements given in Figure 5. For our reference 
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model (Figure 5) the maximum decrease in a v be- 
tween 0 and 10 km was ssl.O MPa. This stress change 
is significant in relation to the stress drops in recent 
earthquakes («2 to 10 MPa) but small in relation to 
the estimated stress magnitude. Therefore the occur- 
rence of an earthquake may be advanced, but the on- 
going stress accumulation would be primarily due to 
tectonic forces. This is consistent with the results of 
numerical simulations of a fault surface with tectonic 
loading by Rydelek and Sacks [1999] which suggest 
that incremental stress changes of several tenths of 
a bar (0.1 MPa) significantly affects the time of an 
earthquake but not its size or its location. 

The location of Mi > 4.0 earthquakes between 
1973 and 1997 was examined to see if there was a 
spatial concentration near regions of glacial retreat. 
In 1979 a large earthquake (St. Elias, Ms = 7.2) 
occurred inland from Icy Bay, and the aftershocks 
associated with this event are evident in Figure 10. 
The source mechanism for this event indicates under- 
thrusting on a northeast dipping plane with a source 
depth of 24 km [Estabrook et a/., 1992]. It is possible 
that the occurrence of this earthquake was advanced 
due to a decrease in a v due to glacial retreat. 

The orientation of postglacial stress release fea- 
tures such as pop-ups and faulting in regions that un- 
derwent extensive late Pleistocene deglaciation sug- 
gest that they were caused by near-surface stresses 
dominanted by radial flexural (fiber or longitudinal) 
stresses near the retreating ice margin [Adams, 1989; 
Adams and Bell , 1991], Determination of focal mech- 
anisms for shallow earthquakes occurring below re- 
cently deglaciated coastal regions could test this hy- 
pothesis. 

7.6. Loading Associated With Offshore 
Sedimentation 

The high uplift rates within the ablation zone of 
the coastal glaciers has led to high erosional rates. 
This then leads to high sedimentation rates in offshore 
basins [Molnia, 1989; Eyles et al, 1991]. Bird [1996] 
speculated that the present localization of crustal 
shortening is a transient effect of Pliocene-Pleistocene 
glaciations which have removed mass from the coastal 
mountains and redistributed it onto the eastern fore- 
arc. He postulated that this may have upset the bal- 
ance between forearc slope and basal traction on the 
plate interface, requiring the onshore part of the fore- 
arc to be shortened and uplifted to restore the dy- 
namic equilibrium. 

Sediment loading would increase cr v and, in a 


thrust earthquake environment, inhibit earthquakes. 
The average thickness of Holocene sediments on the 
south central Alaska continental shelf is about 45 m 
with a maximum thickness of 350 m south of Copper 
River (west of our study region) [ Molnia et al . , 1980]. 
On the basis of vertical stress and pore pressures es- 
timated from the borehole data summarized above, 
this additional loading corresponds to an effective a v 
that ranges from less than 1 MPa to 1.3 MPa. Since 
the end of the Little Ice Age, sedimentation rates near 
Icy Bay have been measured to be more than 1 m/yr 
[Molnia, 1979]. This corresponds to an effective a v of 
almost 1 MPa. Once again these stresses are signifi- 
cant, but we suggest that this would only perturb the 
dominant tectonic stress field. 

8. Summary 

In this study we have estimated crustal deforma- 
tion rates and stress changes in response to ice mass 
changes associated with two major glacial fluctuations 
in the coastal region of southern Alaska between the 
Malaspina and Bering Glaciers. Large uplift and sub- 
sidence on a local scale occurred in response to the ice 
mass redistribution during the Bering Glacier surge. 
An increase in the occurrence of small earthquakes 
beneath the surge reservoir region may be due to a 
decrease in the vertical load in this thrust earthquake 
environment. Glaciological constraints supplemented 
by our geodetic results were used to estimate the ice 
mass redistribution during the surge. Our results sug- 
gest that geodetic data provide new information and 
constraints on ice lowering in the surge reservoir re- 
gion where there is less glaciological information. 

Surface displacements due to the Earth’s viscoelas- 
tic response to retreat of the coastal glaciers dur- 
ing this century were also calculated. The predicted 
yearly uplift rates associated with the rebound pro- 
cesses are generally small. However, in the abla- 
tion zone near Icy Bay the preliminary vertical uplift 
rate, relative to Cape Yakataga, was large enough to 
be measured geodetically, and it may place bounds 
on upper mantle viscosities. The cumulative stress 
changes in the deglaciated region are significant in 
relation to the stress drop in earthquakes but small 
in relation to the estimated tectonic stress levels. Ice 
mass removal could trigger earthquakes earlier than 
would otherwise have been the case, but the primary 
source of stress accumulation is due to ongoing tec- 
tonic processes. 
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Figure 1. Glaciers of the eastern Chugach Mountains, Bering Glacier area, and the western portion of the glaciers 
of the St. Elias Mountains, Yakutat Bay area (modified from Field [1975]). The triangles indicate the location of 
GPS sites (see Table 1). The numbered dots correspond to sites where the ice thinning estimates listed in Table 3 
were made. The Bagley Ice Field occupies a long, narrow east-west trending basin and flows primarily westward 
to the Bering Glacier. Southeast of the Bagley Ice Field is the neve whose principal outlets are the Yahtse and 
Guyot Glaciers in Icy Bay. Farther to the east, a neve in the St. Elias Range is the source region for the Malaspina 
Glacier. The inset shows the region covered by Figure 1. 

Figure 2. A contour plot of the predicted uplift and subsidence (millimeters) of the solid Earth associated with 
the loading/unloading shown in Plate 1. Comparisons of the observed versus predicted displacements are given in 
Table 2. The base map of Plate 1 was used. 


Figure 3. Predicted horizontal elastic displacement rate of the solid Earth associated with recent retreat of 
the coastal glaciers. Disks used to represent the annual unloading rate over the last 30-100 years (circles); the 
magnitude of unloading in meters per year is given from top left to bottom right: -1, -1, -1, -1, -1, -1, -l, -l, -j, _i 
-2, -2, -2, -2, -1, -4, -6, -6, -4, -2, -6, -6, -1 -2, -1, -2, -1, -2, -2, -2, -4, -4, -3, -2, -1.5, -1, -1, -1, -1, -1, -1.5,’ -1.5, -1.5, 
-1.5, -1, -1, -1.5, -1.5, -1.5, -1.5, -1, -1, -1, -1. The primary sites occupied with GPS are given by triangles. The 
Gulf of Alaska coastline is shown by a thick solid line and a simplified outline of the glaciers is given by a thinner 
solid line. 


Figure 4. Contour plot of the predicted elastic uplift and subsidence rate (millimeters per year) of the solid Earth 
associated with the unloading given in Figure 3. 

Figure 5. Finite element grid represention of the subduction zone plate boundary used to calculate the viscoelastic 
response of the Earth to glacial unloading (modified from Cohen [1996]). Given within the key are the assumed 
elastic parameters (E) and viscosity (/j) , in Pa s, for the reference model. 

Figure 6. Schematic of the ice thickness change in the longitudinal profile of a retreating glacier used as input 
to the finite element calculation (after Figure 1 of Schwitter and Raymond [1993]). For the region including the 
peak spanned by a 10 km element in the finite element grid, a maximum value of 300 m was assumed; the other 
elements were scaled on the basis of the profile given. 


Plate 1. Predicted horizontal elastic displacement field (millimeters) of the solid Earth associated with ice transfer 
during the Bering Glacier surge and some thinning due to retreat (black disks, in meters -21, -21) over the 1993-1995 
time frame. The disks are used to represent unloading (red) and loading (blue) in meters (-14, -42, -63, -63, -42, 
-42, -21, -7, -7, -7, -49, -49, -14, -28, -21 -7, -21, -21, -14, -14, -7) and ice loading (14, 14, 14, 14, 14, 14, 14, 14, 14, 
14, 42, 42, 63, 56, 70, 70, 70). The triangles indicate the three sites adjacent to the surge region (DON, ANCX, 
ISLE), one site within «20 km of the surge reservoir region (TIME), and a reference site (VYAK). A satellite 
synthetic aperture radar image of the Bering Glacier and the Bagley Ice Field was used as a base map (see Figure 
1 from Lingle et al. [1993]). The thick solid line near DON and VYAK indicates the approximate position of the 
Gulf of Alaska coast, and the thin solid line indicates the edges of the Bering Glacier and Bagley Ice Field. 

Plate 2. Oblique aerial photograph of the Jeffries Glacier (top) and Bagley Ice Field (bottom), a portion of the 
Bering Glacier, separated by Juniper Island, and a nunatak in the Chugach Mountains, Alaska. The trimline and 
thick lateral moraine draped around the exposed bedrock is indicative of the significant recent thinning that both 
glaciers have experienced. More than 50 m of post-1993 thinning of the Bering Glacier (foreground) has occurred. 
This is due to removal and downglacier transfer of ice during the 1993-1995 surge and postsurge ablation. 



16 

Figure 7. a. Predicted uplift after 1 year (solid line), after 50 years (dash-dot line), and over 100 years (dashed 
line) for the model parameters given in Figure 5. Here the distance from trench refers to the distance from the 2800- 
m contour near the Transition fault [Plafker et ai, 1994], (b.) Predicted horizontal (north-south) displacement 
rate for times as in 7a. A positive displacement rate is toward the trench (« south). 

Figure 8. The Navier- Coulomb failure representation of the stress state of coastal Alaska between Icy Bay and 
Kayak Island. The failure envelope is given for r = 0.4a n and r = 0.85cr„. The Mohr circle representation of 
the effective stresses has been estimated for a depth of 5 km. The estimated stress drop for recent earthquakes is 
compared to the estimated stress with sediment loading over the last 10,000 years, ice mass redistribution during 
the Bering Glacier surge, and retreat this century. 

Figure 9. Earthquakes of M L > 2.5 between 59°N to 61°N and 139°W to 144°W. The dashed rectangles indicate 
the surge reservoir in the Bagley Ice Field and the surge receiving region of the Bering Glacier terminus, (a) 
Earthquakes (79) of Ml > 2.5 from August 1, 1990 through February 28, 1993. (b) Earthquakes (75) of M L > 
2.5 from March 1, 1993, through September 30, 1995. (c) Earthquakes (75) of M L > 2.5 from October 1, 1995, 
through April 30, 1998. The earthquake data are from the National Earthquake Information Center earthquake 
database (// gldss7 . cr . usgs . gov/neis/epic/epic . html) . 

Figure 10. Earthquakes of Ml > 4.0 from 1973 through 1997 between 59°N to 61°N and 139°W to 144°W. The 
largest earthquake to occur during this time period was a M s = 7.2 in February 1979 and the aftershocks associated 
with this event to dominate the seismicity pattern. 



T^ble 1 . Horizontal Site Velocities in a North American Fixed ITRF96 Reference Frame and Vertical Velocities Relative to Cape Yakataga 
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Table 2. Observed Versus Predicted Vertical Displacements(Millimeters), 1993-1995 


Vertical 0 Vertical 6 North 0 North 6 East 0 East 6 

Station Name Predicted Observed Predicted Observed Predicted Observed 


DON 

ANCX 

ISLE 

TIME 

VYAK 


-42.0 

-42.0 

35.4 

32.6 

38.0 

37.8 

20.4 

22.6 

3.7 

0.0 


12.0 

-49.0 

-10.5 

-15.0 

7.4 

8.8 

5.5 

-3.2 

-1.3 

-0.2 


-5.5 

1.8 

1.5 

13.4 

4.1 

11.8 

-4.0 

8.4 

-3.1 

0.0 


^Displacement predicted due to the loading/unloading disk distribution given in Plate 1. 

^ Total displacement estimated from GPS observations made in 1993 and 1995 minus the tectonic strain associated with 
subduction of the Pacific plate beneath interior Alaska. 



Tfeble 3. Constraints on Recent Ice Thinning Indicated by Height Differences Between Glacial Deposits and Nearby Glaciers 

Locatlon Latitude, Longitude, Elevation, m Height Change 4 
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Table 4. Observed Versus Predicted Displacement 
Rates at AMBR relative to Cape Yakataga 


Type 

North 

East 

Vertical 

Observed rate 

6.2 

2.0 

11.7 

Predicted tectonic 

±5 

±5 

±5 

Predicted annual 

0.7 

0.0 

3.1 

5 x 10 20 Pa s 

0.1 

0 

1.4 

5 x 10 19 Pa s 

0.6 

0 

9.3 

5 x 10 18 Pa s 

1.6 

0 

31.2 


Units are millimeter per year. 
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